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We study reaction-diffusion systems where diffusion is by jumps whose sizes are distributed expo- 
nentially. We first study the Fisher-like problem of propagation of a front into an unstable state, as 
typified by the A-l-B — > 2A reaction. We find that the effect of fluctuations is especially pronounced 
at small hopping rates. Fluctuations are treated heuristically via a density cutoff in the reaction 
&J!), rate. We then consider the case of propagating up a reaction rate gradient. The effect of fluctuations 

^ ■ here is pronounced, with the front velocity increasing without limit with increasing bulk particle 

' density. The rate of increase is faster than in the case of a reaction-gradient with nearest-neighbor 

hopping. We derive analytic expressions for the front velocity dependence on bulk particle density. 
\l . Compute simulations are performed to confirm the analytical results. 

'■^ ] Many physical, chemical, and biological systems exhibit fronts which propagate through space. Familiar examples 
^ . range from chemical reaction dynamics such as flames phase transitions such as solidification the spatial 
spread of infections and even the fixation of a beneficial allele in a population "J]. It is thus of great interest to 
I understand the universality classes of fronts which govern what will happen when systems such as these are prepared 
in a spatially heterogeneous manner. These classes determine the selection of propagation speed, the sensitivity to 
particle-number fluctuations, and the stability of the front with respect to deviations from planarity. 

The simplest kind of such a front is that wherein a stable phase replaces a metastable one 0- Here the mean- 
"j^ \ field front velocity is determined via the requirement that there exists a heteroclinic trajectory of the moving-frame 
. steady-state problem (wherein the solution depends only on a; — vt) connecting the metastable phase at -|-oo with the 
I ' stable one at — oo. This type of front is robust with respect to fiuctuations, with power-law corrections in (where 
. A'' is the number of particles per site in the final state) to the mean-field limit A second class is exemplified by 
Ch ' the simple infection model A + B ^ 2 A on a Id lattice (with spacing a) with equal A and B hopping rates Q ; this 
O I process leads in the mean-field limit to a spatially discrete version of the Fisher equation 0] 



^(x) = r<j){x) (1 - (j)ix)) + ^ [<P{x + a) - 2<P{x) + (t){x - a)) . (1) 
>■ ■ 

, Here propagation is into the linearly unstable <f> — Q state, where (j) is the number of A particles at a site. Recent 
■ work IQi, iSi i^J has shown that the front behavior in the stochastic model does approach that of the Fisher equation, 
T— I . where the velocity is selected by the (linear) marginal stability criterion to be 2\frD, albeit with an anomalously 
OO long transient 0(l/i) and anomalously large fluctuation corrections 0(l/ln^A^) . There are also some flndings in 
regard to both front stability in the case of unequal D ^| , and also the scaling properties of front fluctuations jlj] . 
Finally, there are also fronts which have properties intermediate to the previous two classes. 
■»«^ \ In a recent work jl^ Il4j , we introduced a new class of fronts corresponding to propagation into an unstable state 
"j^ ■ up a reaction-rate gradient |l5ltl6j . This type of gradient is present, for example, in systems with an inherent spatial 
inhomogeneity, and also in models of Darwinian evolution |l7l llM Il9ll2(l|. (where the birth rate, which is parallel to 
i' our reaction rate, is proportional to fitness x). We found that the sensitivity to fluctuations in the presence of such 
'■^ a positive reaction-rate gradient is greatly enhanced. In particular, the front velocity diverges with increasing bulk 
^ particle density. As a corollary, the standard reaction-diffusion equation treatment is not useful, as it gives rise to 
O finite-time singularities in the velocity. Also, the velocity is strongly sensitive to details of diffusion, with the increase 
of the velocity with density being qualitatively stronger for a lattice system than in the continuum. 

Given this sensitivity to the precise implementation of diffusion, in this work we turn to the study the effect of 
implementing diffusion via infinite-range hopping, where the size of the jumps is distributed exponentially. Such a 
model has been considered, for example, in the description of the airborne dispersion of seeds, leading to the spread of 
a particular colony of plants. It is also relevant in the evolution context, where the change in fitness due to mutations 
is commonly assumed to be exponentially distributed |2lj |. We will see that even in the absence of a gradient, this form 
of diffusion increases dramatically the effect of fluctuations, at least for small hopping rates. In particular, the naive 
reaction-diffusion formalism predicts a finite velocity in the limit of zero hopping rate, which is clearly unphysical. 
Introducing a reaction-rate gradient again changes the functional dependence of the velocity on density from that of 
the nearest-neighbor hopping studied previously. 

The plan of the paper is as follows. In Section we discuss the gradient-free model, and derive the velocity in 
the limit of infinite density. We show that for fixed hopping rate, the finite density correction formula derived by 
Brunet and Derrida for the Fisher equation is applicable. However, this formula breaks down in the small hopping 
rate limit. We derive an analytical expression for the velocity in this limit. In Section^ we discuss a similar model of 
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Snyder designed to model the spread of colonies, showing that the same physics applies upon the correct mapping of 
parameters. In Section lllll we introduce our reaction-rate gradient model, and after briefly reviewing what is known 
for continuum diffusion and nearest-neighbor hopping, we calculate an analytical approximation to the velocity for 
large density. Finally, in Section llVI we summarize our results and draw some conclusions. 



I. EXPONENTIAL HOPPING FISHER EQUATION 



In this model, the hopping probability has an unbounded range, and decreases exponentially with distance. As in 
the Fisher model 0, the reaction rate is a constant r. In the contimmm limit, the equation describing this model is: 

0(x) = / dse-f^' ((/)(a; + s) + (j){x - s)) - D(i'^cj){x) -f r(j){x) (1 - (j){x)) , (2) 

^ Jo 

and the steady-state equation is: 

I dse-'^' + s) + (t){x - s)) - D0^(j){x) + v(j)' {x) + r0(2;) (1 - ^(x))) = 0. (3) 
2 Jo 

It is useful to convert this equation into a differential equation using the fact that 

0± dse-'^'(j){x±s)= (^(3T^^ dse-^'4>{x ± s) ^ (t>{x). (4) 

Then, acting upon Eq. (j2Jl by O+O- yields 

Z)/3V'+ W + (5) 

As with the standard Fisher equation, this equation has solutions for all velocities, and positive definite solutions for 
all velocities greater than some critical velocity, the so-called marginally stable velocity, vp, which is the asymptotic 
velocity of propagation of all fronts with initial compact support. This can be found from the dispersion relation for 
the leading edge where ~ e"*^^: 

D(3^k^ ^{13^ -k'^){-vk + r)^Q (6) 

The marginally stable velocity, vp, is then given by the requirement that Eq. 10 has a degenerate solution, loading 
to the discriminant condition 

0=4- \DI3'^k'^ + ((3'^ - k^) i-vk -f r)] = 2DI3'^k - (3^v + ivk^ - 2rk. (7) 
dk 



Solving simultaneously Eqs. © and 10 yields, introducing t = ■\/ 0^/3"^ + SD^^r 



= 8^ V r^DP^ 

This has the scaling form vp = f {P^ D / r) where the function f{x) — > 1 for x 00 and f{x) ^ l/{2^/x) as 

a; ^ 0. Thus, for large /3, we recover the usual Fisher answer. What is remarkable is that the velocity has the finite 
limit r//3 as — > 0, so that we have velocity without diffusion! 



A. Calculation of the Velocity for a Small Cutoff 

This anomaly is yet another example of how the reaction diffusion equation, Eq. ^ , provides incorrect information 
about the original stochastic model. A more accurate picture is achieved by studying a cutoff version of the eq uation, 
wherein the reaction is turned off wherever (j){x) is less than some threshold e, of order l/A/'H H El IH lH . This 
captures an essential feature of the original model, namely that the reaction zone always has compact support. Brunct 
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and Derrida ^ have provided a general formula for the correction induced in the velocity due to the cutoff (for small 
cutoffs) for Fisher-like equations. This formula reads 

2(lne)2 

where kp is the degenerate solution of the dispersion relation, Eq. ©. Although derived for second order equations, 
whereas our equation is of third order, nevertheless, as we shall see, it correctly gives the leading order correction for 
the velocity in our case as well. 



(9) 



1. Jump Conditions at the Cutoff Point 



The first task, as for the standard Fisher equation, is to solve the equation for the region beyond the cutoff, where 
4>{x) < e. This will give a set of boundary conditions at the cutoff point, Xc- Due to the third-order nature of our 
equations, and that the derivatives act on the now discontinuous reaction term, these conditions are fairly messly. 
While the solution is continuous at the cutoff point, there is no continuity of the first and second derivatives at this 
point. 

To derive the correct jump conditions, we start from the integral equation, Eq. For x > Xc, the solution is 

(f>x — ee~''''^^~^''\ where kr satisfies the r = dispersion relation 

DP^k^ - {P^ - kl)vkr = 



or, 



Thus, 



dse 



{Xc + s) = 



(3 + kr 



Evaluating the integral equation, Eq. as x — > a;+ gives 



Df3^ 



P + kr Jq 



+ / dse~'^''(j){xc - s) 



D0^e - vkre = 



This, together with Eq. yields 



dse 



(Xc - S) 



P-kr 



Now, let us analyze Eq. Q for x ^ . We get 



n^^) = -'-^^-kre 



Breaking up Eq. © as follows. 



= 



DP^ 



3 r fXa—X 



dse 



2 

- DP^(^ + v(\> + r(\>[\ - 
00^ ' 



2 

- Z)/32(/) -I- -I- r</)(l 
and taking a derivative, we get 



+ s)+ dse~^''(j){x + s) + / dse~f^'(j){x - s) 

J Xc—X 

dssmh{Ps)(l){x + s)+ e^^^^"^-^) 



P + k, 



+ e 



P-kr 



DP^ 



P + kr P - k 



Dp^<j>'{x-) + v(l)"{x-) + r(j)'{x-)(l - 2(l){x-)) = 



or 



-P'^vr + v'^kf + 2vrk'^ + r^kr 2/^^^^ ~ Svrk"^ — Sr'^kr 



krV"^ 



krV^ 



(10) 

(11) 
(12) 

(13) 
(14) 
(15) 



(16) 
(17) 
(18) 
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2. The Modified BD Treatment 

As in the original BD treatment, we divide the range of x < Xc into two regions. In the first region, (j){x) is not small 
compared to 1, but the effect of the cutoff is negligible. In the second region e < 4>{x) << 1. We fix the translation 
invariance by requiring 0(0) = 1/2. Then as e 0, Xc oo. In the first region, we can take the velocity to be vp, 
so that there is a degenerate solution of the dispersion relation. Then, for large x, the dominant solution is 



(t){x) - Axe-'"''= (19) 
In the second region, since the velocity is close to vp, = vp — A, A ^ 1, the general solution is: 

(f>{x) - Be-'''"' sin{hx + C) + Fe''^^^. (20) 



where fc2 < is the third (nondegenerate) root of the dispersion relation and fc^ ^ v A and we can ignore the 
0(A) shift in the real part of k. Matching between the first and the second region requires that B = A/ki, C = 
and i^e^'^'^^^*^^^^'^ 1. Now, in general, we have to enforce three jump conditions, (whose left hand sides are A- 
independent to leading order), with the two free parameters B and F, which is impossible. The only way to make 
things work is to have sm{kiXc) be of the same order as kiCOs{kiXc), in other words kiXc ^ tt — 0(A^/^), which 
is exactly the same condition as in the original BD treatment, where there was one free parameter and two jump 
conditions. Since 



/ 2A 
v"{kp 



(21) 



we immediately recover the BD result quoted above, Eq. 

Examining the BD result, we see that in the limit of 0(3"^ /r 1, kp ^ -^/r/Z? and v"{kp) ^ so that 

v^r^vp- (22) 

In e 



which is of course the Fisher result. On the other hand, when D0^ /r 1, kp ^ P — 0^ ^fDj^r and v"{kp) ~ 
(2r)3/V/3VVi^, and so 

Thus the BD correction diverges as D — > 0. Thus, while for sufficiently small e, the BD correction is correct, for a 
given e, the BD correction fails for small enough D. We show in Figs, ^and^a plot of vp and the BD velocity for 
e = 10~^, compared to the results of an exact numerical calculation. In Fig. [3 it can be seen as predicted that the 
BD treatment does not apply for small D. A calculation in this limit is presented in the next subsection. 



B. Small D, small e limit 

Clearly, in the presence of a cutoff, the velocity should vanish as ^ 0. Let us solve the model in this limit. First, 
let us examine what happens when D = {). Then, for small e we can linearize around the solution 0o = j+^^W^- 
equation reads 

' ^) + "^^^^ " = ^ ^^^^ 

This is equivalent to 

v5' + r5 tanh f — W Ae"'^^ + Be"^ (25) 
V 2u / 



so that 



cosh ^ 



i^^)L5cosh^ r^\\' ^Ae-^^ + Be^^ (26) 
2?; V \2v J J 
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FIG. 1: V vs. large D for exponential Fisher model. Analytical formula for «o ©, and analytical formula for 0, compared 
to numerical results, e = 10~^, /3 = 1. (color online) 




D 



FIG. 2: V vs. small D for exponential Fisher model. Analytical formula for iiq ®, and analytical formula for @, compared 
to numerical results, e = 10~^, (3 ~ 1 (color online) 



and 



2«cosh2(if) 



+B e 



(27) 



A 




V 




B 




V 





We have chosen the limits of integration so that ^(0) = 0, so that the center of the front does not move. What is 
important is the large-x asymptotics of 5: 

|+2e'"/"(-4 + i- I (28) 

-PV V ^ (i) -/3VJ 

We can now use the jump conditions, with — (3 since D = 0, to fix Xc^ A and B. We get, to leading order in e. 



r — Pv 
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B = (29) 
The interesting question is now the behavior at a; — > — oo. The leading asymptotics is 

S - ^c"'^^ I I (30) 

which of course violates the boundary conditions. Thus, there is no solution without D. To leading order in D, we get 
an inhomogeneous term, D(3'^(j)Q, on the left-hand-side of Eq. (|24|l . The inhomogeneous solution, 5d, then satisfies 
the equation 

v6', + r tanh C-^) 5, = H dyG{. - y) ( (31) 
^ \2vJ J-oo V 4i>2 ;cosh3(I|) ^ ^ 



where G is the Green's function for the operator /J^ 



1 



so that 



G{x - y) ^ -e-f'^^-y^ (32) 



dx' cosh^ {^] I dyG{.' - y) ( ^^M, (33) 



i;cosh2(i|) 7o V2V7-00 " ' 4«2 ;cosh3(IH) 

We need the asymptotic behavior of 5d for large x. For r/v ^ P, the integral is dominated by the region of x' large, 
y 0. Thus, 



4g-rx/. rx rx'/v rco ^-f3x' / ^2 2 X / ^^2^2 X giuh (^) 

' dx' — ^ — / dy—— 1 + f3y+ + ' ' ' 



V 



2(i \ ' " 2 ■ ■ 7 V 4t;2 ; cosh^ (§f ) 



2v -B V 8r2 

We now have to again solve the jump conditions with this new contribution. The coefficient A above is now modified 
and includes a term which, up to linear order in Pv/r, reads 



2 \ r ^ 

The condition for a solution is that this cancels the A we found above, so that 



1 + ^ ^ ew/3e'^^= (36) 
z \ r J 

or 

n - 2et; i3v/r\n{ernT-pv)) _ 2t;r i-p^/r ( _!__\ ^ ^ /07^ 



Thus, for very small D, the velocity is equal to DP/{2e), which is reminiscent of the behavior of evolution models for 
very small mutation rates [l^ . The comparison between our analytic approximation and an exact numerical solution 
is shown in Fig. 
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FIG. 3: Comparison of our analytic approximation, Eq. 13711 and an exact numerical solution of Eq. Parameters are 

e — 10"^, (i = r — 1. (color online) 



II. THE SNYDER DISCRETE-TIME MODEL 



Recently, Snyder 7] introduced a model of colony spreading which, in one variant, involved an exponentially- 
distributed hopping similar to the model defined above. The essential difference between her model and ours is that 
hers was a discrete-time model. In each time step, all the offspring performed a hop and the parental generation was 
removed. The number of offspring at a given site was given by a local logistic growth law, similar to that incorporated 
in the Fisher model. Snyder performed numerical simulations and measured the velocity of propagation, both for 
the stochastic model, and for the corresponding (uncutoff) reaction-diffusion system, and found a difference between 
these two velocities. Due to its close correspondence to the present model under investigation, it is useful to derive 
analytically the uncutoff velocity and the BD approximation to the cutoff velocity, so as to make clear the mapping 
between the Snyder model and ours. 

As always, to derive the uncutoff "Fisher" (marginally-stable) velocity, it is enough to consider the linearized version 
of the Snyder model, which reads 



r(3 



-p\y-x\ 



dy 



(38) 



where rg is the average number of offspring per individual and 0* (y) is the number of individuals at site y at (integer) 
time t. In Fourier space our equation reads: 



which, starting from a (5-function initial condition gives 

4>\k) = 

or, Fourier transforming back, 

0*(x) 



/32 +/C2 



rsf3^ 







J —oo 


-|-fc2_ 



,dk 
2^ 



(39) 



(40) 



(41) 



We want to calculate the velocity, so we are interested in (j){x — vt), where we have to choose v such that this is 
independent of t for large t. This gives us a saddle-point integral 



°° tin SS£_ 



dk 

27r 



The saddle point is at fc, where 



(42) 



(43) 
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The dominant contribution to the integral is given by evaluating the integrand at the saddle, giving 

expftflnf-^V^^*-)) (44) 



If this is to be independent of t, the term in the exponential must vanish: 

1- f#4) + A - (45) 
Clearly, k is pure imaginary and proportional to /3, so we write 

fc, = ial3 (46) 



so that a depends only on r and satisfies 



and the velocity is 



It is reassuring that this formula reproduces the velocity measured by Snyder for the one set of parameters presented 
in her paper. For rs near 1, w ~ 2y/ (r — l)/2//3, while for large rs, v ^ lii(rs)/(3. Of course, on dimensional grounds 
this is reasonable, since w is a velocity per round, which has units of length, and rs is dimensionless. We see that 
rs near 1 corresponds to the Fisher limit, equivalent to the large D0^ /r limit of our model, since the growth rate of 
the population is — 1, so that small — 1 corresponds to a large value of our dimensionless control parameter. 
On the other hand, for large rg, the models differ since the diffusion never goes away entirely in the everyone hops 
Snyder model. It is also interesting to note that in the Snyder model, the only effect of (3 is to set the velocity scale, 
as opposed to the more complicated role of [3 in our model. 

In fact, there is another way to solve equation (|38|) . We assume that the dependence of in i and y is 

cl>\y)^cj>{y-vt), (49) 

and 

^*(y) =e-"(^-*''). (50) 

putting IHOI) in ijSHll yields: 



/32 - a2 ■ 

Taking the derivative by a of H51|l (according to the marginal stability criterion), and dividing it by H51f) yields: 



(51) 



lap 

VF = (52) 



Eqs. 1)51(1 and ((52|l are seen to be equivalent to Eqs. 147|l and H48|l upon defining a = ap/(3. 
We can eliminate a_F to obtain a direct relationship between r and fF as follows: 



-1 + + 1 

ap = — • (53) 

Vp 



so 



2(^/1 + 4/32 ~l)eV4/3^+i-i 



(54) 



v'pf3^ 

II! \ 1 1 

The advantage of this second approach is we can immediately write down the BD correction, = vp — '" "^^^^a ""^ ■ 
A graph of Snyder velocity with and without the correction is shown in Fig. 01 We see that the larger rs is, the larger 
the correction according to BD is, since as discussed above, increasing rg corresponds to decreasing the strength of 
diffusion in our model. 
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FIG. 4: Snyder model: vp vs. rg according to 1521 1 and the BD correction vs rs- /3 — 0.5 (color online) 



III. EXPONENTIALLY DISTRIBUTED HOPPING WITH A REACTION-RATE GRADIENT 

In a previous work [isl , we studied the case of frontspropagating into an unstable state up a reaction-rate 
gradient. We focused again on the A + B ^ 2 A reaction with no A particles and an initial mean number TV 
of B particles at all sites past some initial xq, but with a reaction probability that depended linearly on spatial 
position. This type of gradient would be a natural consequence of spatial inhomogeneity, or could be imposed via a 
temperatu re g radient in a chemical reaction analog. Also, this type of system arises naturally in models of Darwinian 
evolution |l7l ll9j, (where fitness x is the independent variable; the birth-rate, akin to the reaction-rate here, is 
proportional to fitness). The naive equation describing such a model is the Fisher equation with a reaction 
strength r — ra{x) varying linearly in space 

ra{x) = max(rnii„, tq + ax) . (55) 

where r^i„ is introduced to insure that the reaction rate stays positive far behind the front, and has no effect on 
the velocity. This model gives rise to an accelerating front. We also introduced a quasi-static version of the model, 
wherein the reaction rate function moves along with the front: 

rq{x) = max(rinm, fo + a{x - Xf)) , (56) 

with Xf is the instantaneous front position. This quasi-static problem should lead to a translation- invariant front 
with fixed speed Vq(fQ,a). Although important on its own, one might also try to view the quasi-static problem as 
a zeroth-order approximation to the original model, (the absolute gradient case) , where by ignoring the acceleration, 
one obtains an adiabatic approximation to the velocity v(t;ro,a) ~ Vq{fo{t), a) with fo(t) = tq + axf{t). In both 
models, fluctuations become crucial due to the reaction gradient and the presence of the gradient leads to a new 
class of fronts. One characteristic of this new class is the divergence of the front velocity with N. We found, that to 
leading order, the velocity of the front in the continuum limit diverges as ln.^^^{N), and to leading order on a lattice, 
the velocity diverges as vlniV. It should be noted that in both cases the leading order does not yield an accurate 
solution, and the next order correction must be taken into account. 

Given that the nature of the divergence of the velocity with N depends on the microscopic implementation of diffu- 
sion (continuum versus lattice) , it is natural to investigate this question for our model with exponentially distributed 
hopping. Here, we chose to work on a lattice (with spacing a); we will see in the end that the results here are not 
sensitive to the presence of the lattice. The model we study is: 



dt a^e'r{e'y + 1) 



^(e-^^"(0,+,+0._,))-2 
i=i 



- 1 



+ r{i)Ml~<f>,)0{p,-e), (57) 

where 7 == /3a is the rate of exponential falloff of the hopping between successive lattice sites. It is easy to verify 
that this model reproduces continuum diffusion with coefhcicnt D for sufficiently smooth fields (pi . We choose to focus 
on the quasi-static problem, as the presence of a steady-state solution makes the problem analytically tractable. The 
steady-state solution on the lattice has the Slepyan [2J| form 

(j),{t) = (j){t-ia/v) (58) 
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so that each lattice point experiences the same history, with a time shift. 
—v(t — ia/v), in terms of which 



We define the continuous variable 



= D 



(e-^ - 1)^ 



a2eT (eT + 1) 
+ r(z)0(z)(l - ^{z))9{p, -e) + v(j)'{z), 



{(j){z + aj) + (j){z - aj))] 



(^) 



(59) 



We wish to solve this equation for small e, assuming that v will be large in this limit. Relying on our previous 
analysis of the nearest-neighbor hopping problem, we expect that the leading order solution for the velocity comes 
from the region of the front where (j) is small, so the nonlinear (j)^ term can be dropped. We assume 13., 2^ a 
WKB-type solution 



e^" , and expand Sz±aj into a Taylor series, so equation H59() becomes: 
{e' - 1)3_ 

T) 

2 



= D 



a?e^{ei + 1) 
1 



1 



+ az + vS' . 



1 



el - 1 



After some algebra we get from (|6U|) : 

AD 

= 



sinh2(aS"/2) 



=27 



1 - 26-)- cosh(S''a) 



vS'. 



(60) 



(61) 



As in the nearest-neighbor hopping problem, the only way to match to the post-cutoff solution is to require that the 
front be close to the classical turning point. In order to find the turning point, we need to equate the derivative of 
(|60|1 with respect to S" to zero. Doing so we get: 



: 



2D 



sinh(5^a) 



(e27 + 1 - 26-^ cosh(5»)2 



(62) 



where S'^ is the value of S' at the turning point. For large 7, Eq. (|62|l indeed matches the nearest-neighbor hopping 
result. For large v, the denominator of the first term in Eq. H62|l has to be close to in order to balance the second 
term. As the denominator vanishes if aS'^ — —7, this gives 



-1 
a 



2a^ 



_27 (£1^1)1 
sinh(7) 



(63) 



7/a, SI is only weakly dependent on 
5*^ = -1 y^D/^ for a ^ and 



1 



which is correct for v ^ Dj/a. From (|63|l one can obtain that, for fixed /3 
a for < a < 1 as long as f3 is not too large, /? < 2. For example, for f3 
— 1 + LOOIOyZ-D/w for a ~ 1. This is reasonable, since the long-range nature of the hopping (for not too large /?'s), 
smooths over the lattice structure. 

The fact that \Sl\ is bounded by (3 is the unique feature of our exponentially-distributed hopping. We remind the 
reader that for standard continuum diffusion, ~ —v/{2D), and so \S'J grows unboundedly with v, while for nearest- 
neighbor hopping, \Sl\, though not linearly dependent, still grows logarithmically with v. The faster the growth of 
\S'^\, the weaker the dependence of the velocity on ln(e). This confirms our initial intuition that the exponentially 
distributed hopping model should be more sensitive to fiuctuations that even the nearest-neighbor hopping model. It 
also reiterates why the lattice parameter a is not important (for /? not large) , since the rate of exponential falloff of 
(j) is bounded by (3, and so never gets too large as to be affected by the lattice. 



Since the turning point is close to the cutoff point, the dominant contribution to the value of i 
the value of S at the turning point. We now want to find S^. This is given as : 



is e^' , where S^, 



S ^ 



dzS' 



s' 



si I 

dS'S'i — ) 
a 



2D 



{e'^ - l)4sinh(7) 



a (e2T - 2e'rcosh(S''a))2 
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D 



D 



1 



{e'' - 1)" 
(e^ - if 



S 



e27 + 1 _ 2eT cosh(S'» 



1 



a(e27 - 1) 



In 



ef - e 



aSL 



2a 



(64) 



To leading order, 
point, 



(zc) = e'^*. In order to get the correction for S*,, we write, in the vicinity of the turning 



(z) = e^-^V(^) 



(65) 



Equation H65(l smooths the variation between lattice points in the vicinity of the turning point, so we can expand 
ipi^z) in a Taylor series. This gives 



= £>■ 



2eT(e^ + 1) 



^(e-'^J'(e^-'^J'V(z + aj) + e-^'-^^^iz - aj))) 



2D- 



(e-y - 1)3 
a2eT(eT + 1) 



eT - 1 



+ (ro + q;z)-0(z) + vipiz)S' + vip'iz) 



DiP{z] 



(e^ - 1)3 



Dil){z) 



a?e'^{e^ + 1) 
(e-^ - 1)' 



1 



1 



a?e^{e'i + 1) 



Di)'{z) 



026^(6-^ + 1) 

(e^ - 1)3 



+ D^l)"{z) 



27+Si _ 1 e7 _ 1 



a2e'^(eT + 1) 
+ (rg + az)?/;(z) + v'ip{z)S' + vijj'{z). 
After some algebra, we get 

(e'^ - 1)-* [(e^-^ + 1) cosh(S'» + 2e'' cosh^(S'» - 4e'^] „ 



(66) 



= D- 



(e27 + 1 - 26-/ cosh(5'»)= 



riz) 



+ a{z — z^,)ip{z) 
This is the Airy equation. The solution of ()67|l is 

^p{z) = Ai {z„ — z) 



a 



(e2T + l-2e'^ cosh(S'»)3 



D(ert - 1)4 [(e27 + 1) cosh(S':a) + cosh2(S':a) - \e^\ 



(67) 



(68) 



This gives us the distance from the turning point to the zero of 0. The first zero of the Airy hmction is at —2.338, so 
that the distance, I, is 



i = 2.338 



a(e2T + 1 - 2eT cosh(5»)3 



p{ei - 1)4 [(e27 + l)cosh(S'» + 2e7 cosh2(S':a) - ^ei\ 
This gives us an addition contribution to S*, of S'^t. Adding this to Eq. (|64() yields: 



(69) 



In 



D 



{e'< - 1)4 
(e-'- 1)4 



6^7 + 1 - 2e'>' cosh(S'» 



In 



o7+aSt _ 1 



a(e27 - 1) I eT 
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FIG. 5: V vs. ln{N) for D = 1, a = 1, a = 1, q = 0.1, B = 1. Numerical simulations are compared to the 1st order formula, 
Eq. 17U1 and the 2nd order formula, Eq. I(i4ll . (color online) 



+ 7^(^1)'« - 2.33851 f—^— 

2a \D{e'i - 1)^ 



(6^-^ + 1-26^ cosh(5»)3 

[(e27 + 1) cosh(S'» + 2e')' cosh^ {S'^a) - Ae^] 



(70) 



Again, this solution matches our previous solution for /3 >> 1 [l^ Il4|. In the continuum hmit, which as we noted 
above is accurate for /3 < 2, this equation becomes 



S ^ 



- (Si)' 



DI3\ f(3 + S: 
■ In 



2a 



— (5l)\- 2.33851 
2a 



/52 + 3(502 



(71) 



Substituting the continuum limit of our expression for 51 in the above yields, and expanding for large v yields 



e 2a 



2.338 



2^/^VP 



which is still fairly messy. To test these formulas, we present in Fig. [31the velocity versus ln(l/e), comparing between 
Eqs. H7U|) and H64I) and numerical results from direct integration of the time-dependent equation. We see that the 
agreement between theory and simulation is quite good, and that the correction term is not negligible for this range 
of ln(l/e). 

The first interesting thing to note about our analytic result is that asymptotically, for small cutoff, the velocity is 
proportional to ln(l/e), with a coefficient independent of D. This is reminiscent of the "velocity without diffusion" 
we saw in the zero-gradient ca se in th e absence of a cutoff. We can see this point clearly in Fig where we graph 
?;/ln(l/e) as a function of 1/ ^/hiJTJe) for Z? = 1 and D = A. It is clear that the two graphs are converging to the 
same value of /3^/(2a) = 0.2. 

In Fig. [7|we present results for v versus a, again comparing the analytic formulas Eqs. H70|) and (|64|l to the results 
of direct simulation. Again the agreement is very satisfying. For large a the "correction" term is dominant and the 
velocity grows as. 



0.1452 



ln^{l/e)a^/^D^^^ 



(73) 



Unfortunately, this asymptotic result is only valid for extremely large a 3> ln^(l/e). 

The dependence of the velocity on the diffusion constant D is presented in Fig. [SJ where we again present a compar- 
ison with our theoretical prediction. It is seen that as D grows, v/ D decreases, and so our analytic approximation for 
51 becomes increasing less reliable. Further analysis shows that in fact for very large D, 51 « v/{2D) ^ 1, and the 



13 



0.8 r 



0.6 



o Numerical D=l • 
• Numerical D = 4 



C 0.4h 

CEOO O 

0.2 - 



In 



0.2 



FIG. 6: v/ln{N) vs. ln{N) for D = 1, a = 1, a = 1, a = 0.1, ro = 1. The data presented are from numerical simulations of 
the cutoff deterministic equation, (color online) 
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FIG. 7: V vs. a for D = 1, (3 — 1, a = 1, ro = 1, ln{N) = 25. Numerical simulations are compared to the 1st order formula, 
Eq. (1701 and the 2nd order formula, Eq. I64II . (color online) 



calculation reverts to that of the standard continuum diffusion presented in Ref. where v ~ /(a, e)D^/^ , and the 
prefactor / ~ (24aln(l/e))^/'^ for e — > 0. We can verify this result by replotting the data in Fig. |51 this time showing 
u/D^/^, which is seen to be consistent with an approach to a constant close to (24aln(le))^/'^ — 3.91. This reversion 
to continuum diffusion for large D is reasonable, since if diffusion is fast enough, it is irrelevant how it is implemented. 
For extremely small D our calculation becomes unreliable, since there one is not allowed to truncate to an Airy 
equation. We expect, similar to what we occurs in the evolution problem, that the velocity will be proportional to 
D/e in this limit. 

Lastly, In Fig. EH we can see a comparison between (|7()ll . H64|l and numerical results for v vs. the rate of falloff 
of the hopping distribution, [3. For large (3, the problem reverts to the nearest neighbor hopping model, so v should 
approach a constant in that limit, consistent with the data presented. For small 7, again the "correction" term is 
dominant and we recover the large a result, Eq. H73|l . with v diverging as 1//3. We therefore plot v(3 versus /3 in Fig. 
111! where we see that the data is consistent with vj3 approaching the constant 0.14521n^(l/e)a^/'^i?"'^/^ = 19.55 for 
small p. 

The last task before us is to test our cutoff theory is a good approximation for the stochastic case. The analytical 
procedure done above is referring to the case for which the front position Xf is defined to by 4>{xf) = 1/2. For the 
stochastic case this procedure is ill-defined, since </> fluctuates. Rather, we choose to define the front by 

k 

Rather than redo the theory for this definition of the front, we chose the expedient of comparing the the stochastic 
results to numerical results that also define the front position as the sum of 4>k^ which amounts to a shift in tq. The 
comparison is shown in Fig. 1121 

As a closing remark, we note that one of the most interesting aspects of the above calculation (and the previously 
published calculations for the nearest neighbor hopping model) is that the result does not at all depend on form of 
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FIG. 8: V vs. D for P = 1, a = 1, a = 0.1, ro = 1, ln(l/e) = 25. Numerical simulations are compared to the 1st order formula, 
Eq. 1701 and the 2nd order formula, Eq. I(i4ll . (color online) 
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FIG. 9: v/D'^/^ vs. D ioT fi = 1, a = 1, a = 0.1, ro = 1, ln{l/e) — 25. Data presented are from numerical simulations of the 
deterministic cutoff equation, (color online) 



the solution past the cutoff point; the mere existence of a cutoff is enough to force the system to the WKB turning 
point and hence fix the velocity. 




FIG. 10: V vs. /3 for D — 1, a = 1, a — 0.1, ro — 1, ln{l/e) — 25. Numerical simulations are compared to the 1st order 
formula, Eq. II7UII and the 2nd order formula, Eq. 1641 . (color online) 
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FIG. 11: w/3 vs. P for D = 1, a — 1. a = 0.1, ro — 1, ln{l/e) — 25. Data presented are from numerical simulations of the 
deterministic cutoff equations, (color online) 
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FIG. 12: Comparison between stochastic results and numerical results for v vs. ln{N) for D = 1, a = 1, a = 0.1, ro = 1. 
Xf ~'}2 4'k- (color online) 



IV. SUMMARY 



We have investigated herein reaction-diffusion systems in which the hopping probabihty exponentially decays vifith 
distance, focussing on the fluctuation induced anomalies seen in the same systems with continuous diffusion and 
nearest neighbor hopping. As in these previously studied cases, we probe the sensitivity to fluctuations by studying 
the dependence of the steady-state velocity on a cutoff in the reaction term when the density drops below a cutoff 
of the order of one particle per site. We flrst studied this model with no gradient, showing that, in the absence of a 
cutoff, the velocity does not vanish for small D. We showed that the BD correction for velocity due to the presence 
of a cutoff diverges in the case of small D, and calculate that the velocity actually vanishes linearly for small D 
in the presence of cutoff. Our model is similar to a discrete-time model describing the spread of colonies, and we 
show the same generic features apply to this model as well. We then studied the effect of introducing a quasi-static 
gradient into our model. Here, even for continuum diffusion and nearest-neighbor hoppings, fluctuation effects lead 
to a divergence of the velocity with increasing particle density TV. We found that this phenomenon is enhanced by 
the exponential distributed hopping, so that the velocity diverges more strongly, as ln(l/e). In fact, for long-range 
hopping, /? <C 1, the velocity is proportional to ln^(l/e). Our analytical work was confirmed by direct simulation of 
the cutoff deterministic equation, as well as by comparison to the original stochastic model. 
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APPENDIX: NUMERICAL SIMULATIONS 

In the body of the paper, we have presented results from dhect numerical simulations of both the deterministic 
cutoff reaction-diffusion equation and the stochastic particle model. Here we briefly present some relevant details of 
the simulation methods, especially in reference to treating in an efficient manner the long-range nature of the hopping. 

1. Deterministic Equation 

The simulations are essentially standard, using an Euler method time step. The only subtlety is in handling the 
hopping term efhciently. A naive treatment would involve calculating the transfer of density from every pair of sites, 
which is a prohibitively expensive 0{L^) operation, where L is the spatial extent of the lattice. 

To solve this difficulty, consider the density transfered to site i from all the sites to the left; i.e., 1, 2 ... i — 1. This 
transfered density, which we denote Li is given by 

i-l 

L,=^0,e-'^(*-^") (A.l) 

Li satisfies a simple recursion relation: 

L,; = -t- 4-i)e-^ (A.2) 

Thus, in one pass we can calculate how much density is transferred to every site from all its left neighbors. The 
density transferred from the right neighbors is done similarly, using 

L 

and the recursion relation 

i?, +0,+i)e-T (A.4) 

and making a leftward pass over the sites. The problem is thus reduced to an 0{L) problem. 

2. Stochastic Simulation 

Our basic technique for simulating the stochastic model is to treat all the particles on a given site in " bulk" 
The number of particles that participate in any given process (birth, death and hopping) is given by a binomial 
distribution, and so can be determined by drawing a binomial deviate. The simulation performs in parallel first a 
hopping step, followed by a reaction step. In the reaction step, the number of B particles which transform into A's 
at site X is again a binomial deviate, drawn from B{Nb{x), 1 — (1 — r(x)(ii)^'*^^''). Replacing the distribution by its 
expected value, and setting Nb{x) = N — Na{x), and defining = Na/N gives Eq. A dt small enough so that 
less than 10% of the A, i?'s at a site hop and/or react in one time step is sufficient; smaller values do not alter the 
results. 

Again, hopping in our model provides a challenge, since we cannot afford to draw a binomial deviate for every pair 
of sites. Rather, every time step we first determine the number of particles leaving that site due to the hopping, by 
drawing a single binomial deviate. We then determine how many of these move to the right, by drawing a second 
deviate. Of those moving to the left (right), we determine how many move to the nearest neighbor, by drawing a third 
deviate, and remove this number from the pool of left (right) movers. Then, if any particles remain in the pool, we 
determine how many move to the second nearest neighbor, removing these from the pool, continuing in this manner 
till the pool is exhausted. The number of deviates we need to choose is thus fixed (on average) by 7, independent of 
L. 
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